The present invention relates to processing of seismic data representative of subsurface features in the earth and, more particularly, to improved methods of processing seismic data using limited Radon transformations to remove unwanted noise from meaningful reflection signals.
Seismic surveys are one of the most important techniques for discovering the presence of oil and gas deposits. If the data is properly processed and interpreted, a seismic survey can give geologists a picture of subsurface geological features, so that they may better identify those features capable of holding oil and gas. Drilling is extremely expensive, and ever more so as easily tapped reservoirs are exhausted and new reservoirs are harder to reach. Having an accurate picture of an area""s subsurface features can increase the odds of hitting an economically recoverable reserve and decrease the odds of wasting money and effort on a nonproductive well.
The principle behind seismology is deceptively simple. As seismic waves travel through the earth, portions of that energy are reflected back to the surface as the energy waves traverse different geological layers. Those seismic echoes or reflections give valuable information about the depth and arrangement of the formations, some of which hopefully contain oil or gas deposits.
A seismic survey is conducted by deploying an array of energy sources and an array of sensors or receivers in an area of interest. Typically, dynamite charges are used as sources for land surveys, and air guns are used for marine surveys. The sources are discharged in a predetermined sequence, sending seismic energy waves into the earth. The reflections from those energy waves or xe2x80x9csignalsxe2x80x9d then are detected by the array of sensors. Each sensor records the amplitude of incoming signals over time at that particular location. Since the physical location of the sources and receivers is known, the time it takes for a reflection signal to travel from a source to a sensor is directly related to the depth of the formation that caused the reflection. Thus, the amplitude data from the array of sensors can be analyzed to determine the size and location of potential deposits.
This analysis typically starts by organizing the data from the array of sensors into common geometry gathers. That is, data from a number of sensors that share a common geometry are analyzed together. A gather will provide information about a particular spot or profile in the area being surveyed. Ultimately, the data will be organized into many different gathers and processed before the analysis is completed and the entire survey area mapped.
The types of gathers typically used include: common midpoint, where the sensors and their respective sources share a common midpoint; common source, where all the sensors share a common source; common offset, where all the sensors and their respective sources have the same separation or xe2x80x9coffsetxe2x80x9d; and common receiver, where a number of sources share a common receiver. Common midpoint gathers are the most common gather today because they allow the measurement of a single point on a reflective subsurface feature from multiple source-receiver pairs, thus increasing the accuracy of the depth calculated for that feature.
The data in a gather is typically recorded or first assembled in the offset-time domain. That is, the amplitude data recorded at each of the receivers in the gather is assembled or displayed together as a function of offset, i.e., the distance of the receiver from a reference point, and as a function of time. The time required for a given signal to reach and be detected by successive receivers is a function of its velocity and the distance traveled. Those functions are referred to as kinematic travel time trajectories. Thus, at least in theory, when the gathered data is displayed in the offset-time domain, or xe2x80x9cX-Txe2x80x9d domain, the amplitude peaks corresponding to reflection signals detected at the gathered sensors should align into patterns that mirror the kinematic travel time trajectories. It is from those trajectories that one ultimately may determine an estimate of the depths at which formations exist.
A number of factors, however, make the practice of seismology and, especially, the interpretation of seismic data much more complicated than its basic principles. First, the reflected signals that indicate the presence of geological strata typically are mixed with a variety of noise.
The most meaningful signals arc the so-called primary reflection signals, those signals that travel down to the reflective surface and then back up to a receiver. When a source is discharged, however, a portion of the signal travels directly to receivers without reflecting off of any subsurface features. In addition, a signal may bounce off of a subsurface feature, bounce off the surface, and then bounce off the same or another subsurface feature, one or more times, creating so-called multiple reflection signals. Other portions of the signal turn into noise as part of ground roll, refractions, and unresolvable scattered events. Some noise, both random and coherent, is generated by natural and man-made events outside the control of the survey.
All of this noise is occurring simultaneously with the reflection signals that indicate subsurface features. Thus, the noise and reflection signals tend to overlap when the survey data is displayed in X-T space. The overlap can mask primary reflection signals and make it difficult or impossible to identify patterns in the display upon which inferences about subsurface geology may be drawn. Accordingly, various mathematical methods have been developed to process seismic data in such a way that noise is separated from primary reflection signals.
Many such methods seek to achieve a separation of signal and noise by transforming the data from the X-T domain to other domains. In other domains, such as the frequency-wavenumber (F-K) domain or the time-slowness (tau-P), there is less overlap between the signal and noise data. Once the data is transformed, various mathematical filters are applied to the transformed data to eliminate as much of the noise as possible and, thereby, to enhance the primary reflection signals. The data then is inverse transformed back into the offset-time domain for interpretation or further processing.
For example, so-called Radon filters are commonly used to attenuate or remove multiple reflection signals. Such methods rely on Radon transformation equations to transform data from the offset-time (X-T) domain to the time-slowness (tau-P) where it can be filtered. More specifically, the X-T data is transformed along kinematic travel time trajectories having constant velocities and slownesses, where slowness p is defined as reciprocal velocity (or p=1/v).
Such prior art Radon methods, however, typically first process the data to compensate for the increase in travel time as sensors are further removed from the source. This step is referred to as normal moveout or xe2x80x9cNMOxe2x80x9d correction. It is designed to eliminate the differences in time that exist between the primary reflection signals recorded at close-in receivers, i.e., at near offsets, and those recorded at remote receivers, i.e., at far offsets. Primary signals, after NMO correction, generally will transform into the tau-P domain at or near zero slowness. Thus, a mute filter may be defined and applied in the tau-P domain. The filter mutes, i.e., excludes all data, including the transformed primary signals, below a defined slowness value Pmute.
The data that remains after applying the mute filter contains a substantial portion of the signals corresponding to multiple reflections. That unmuted data is then transformed back into offset-time space and is subtracted from the original data in the gather. The subtraction process removes the multiple reflection signals from the data gather, leaving the primary reflection signals more readily apparent and easier to interpret.
It will be appreciated, however, that in such prior art Radon filters, noise and multiple reflection signals recorded by receivers close to the gather reference point (xe2x80x9cnear trace multiplesxe2x80x9d) are not as effectively separated from the primary reflections. The lack of separation in large part is an artifact of the NMO correction performed prior to transformation. Because NMO correction tends to eliminate offset time differences in primary signals, primary signals and near trace multiple signals both transform at or near zero slowness in the tau-P domain. When the mute filter is applied, therefore, the near trace multiples are muted along with the primary signal. Since they arc muted, near trace multiples are not subtracted from and remain in the original data gather as noise that can mask primary reflection signals.
Radon filters have been developed for use in connection with common source, common receiver, common offset, and common midpoint gathers. They include those based on linear slant-stack, parabolic, and hyperbolic kinematic travel time trajectories. The general case forward transformation equation used in Radon filtration processes, R(p,xcfx84)[d(x,t)], is set forth below:       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢                        ∫                      -            ∞                    ∞                ⁢                  xe2x80x83                ⁢                              ⅆ            t                    ⁢                      xe2x80x83                    ⁢                      ⅆ                          (                              x                ,                t                            )                                ⁢                      δ            ⁡                          [                              f                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                ⁢                      xe2x80x83                    ⁢                      (                          forward              ⁢                              xe2x80x83                            ⁢              transformation                        )                              
where
u(p,xcfx84)=transform coefficient at slowness p and zero-offset time
xcfx84d(x,t)=measured seismogram at offset x and two-way time t
p=slowness
t=two-way travel time
xcfx84=two-way intercept time at p=0
x=offset
xcex4=Dirac delta function
f(t,x,xcfx84,p)=forward transform function
The forward transform function for linear slant stack kinematic travel time trajectories is as follows:
f(t,x,xcfx84,p)=txe2x88x92xcfx84xe2x88x92px
where                               δ          ⁡                      [                          f              ⁡                              (                                  t                  ,                  x                  ,                  τ                  ,                  p                                )                                      ]                          =                  xe2x80x83                ⁢                  δ          ⁡                      (                          t              -              τ              -              px                        )                                                            =                      xe2x80x83                    ⁢          1                ,                  xe2x80x83                ⁢                              when            ⁢                          xe2x80x83                        ⁢            t                    =                      τ            +            px                          ,        and                                          =                      xe2x80x83                    ⁢          0                ,                  elsewhere          .                    
Thus, the forward linear slant stack Radon transformation equation becomes       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢              ⅆ                  (                      x            ,                          τ              +              px                                )                    
The forward transform function for parabolic kinematic trajectories is as follows:
f(t,x,xcfx84,p)=txe2x88x92xcfx84xe2x88x92px2
where                               δ          ⁡                      [                          f              ⁡                              (                                  t                  ,                  x                  ,                  τ                  ,                  p                                )                                      ]                          =                  xe2x80x83                ⁢                  δ          ⁡                      (                          t              -              τ              -                              px                2                                      )                                                            =                      xe2x80x83                    ⁢          1                ,                  xe2x80x83                ⁢                              when            ⁢                          xe2x80x83                        ⁢            t                    =                      τ            +                          px              2                                      ,        and                                          =                      xe2x80x83                    ⁢          0                ,                  elsewhere          .                    
Thus, the forward parabolic Radon transformation equation becomes       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢              ⅆ                  (                      x            ,                          τ              +                              px                2                                              )                    
The forward transform function for hyperbolic kinematic travel time trajectories is as follows:
f(t,x,xcfx84,p)=txe2x88x92{square root over (xcfx842+p2x2)}
where                               δ          ⁡                      [                          f              ⁡                              (                                  t                  ,                  x                  ,                  τ                  ,                  p                                )                                      ]                          =                  xe2x80x83                ⁢                  δ          ⁡                      (                          t              -                                                                    τ                    2                                    +                                                            p                      2                                        ⁢                                          x                      2                                                                                            )                                                            =                      xe2x80x83                    ⁢          1                ,                  xe2x80x83                ⁢                              when            ⁢                          xe2x80x83                        ⁢            t                    =                                                    τ                2                            +                                                p                  2                                ⁢                                  x                  2                                                                    ,        and                                          =                      xe2x80x83                    ⁢          0                ,                  elsewhere          .                    
Thus, the forward hyperbolic Radon transformation equation becomes       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢              xe2x80x83            ⁢              ⅆ                  (                      x            ,                                                            τ                  2                                +                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
A general case inverse Radon transformation equation is set forth below:       ⅆ          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        p            ⁢                        ∫                      -            ∞                    ∞                ⁢                  xe2x80x83                ⁢                              ⅆ            τ                    ⁢                      xe2x80x83                    ⁢                      ⅆ                                          ρ                ⁡                                  (                  τ                  )                                            *                                ⁢                      u            ⁡                          (                              p                ,                τ                            )                                ⁢                      δ            ⁡                          [                              g                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                          
where
g(t,x,xcfx84,p)=inverse transform function, and
xcfx81(xcfx84)*=convolution of rho filter.
The inverse transform function for linear slant stack kinematic trajectories is as follows:
g(t,x,xcfx84,p)=xcfx84xe2x88x92t+px
Thus, when xcfx84=txe2x88x92px, the inverse linear slant stack Radon transformation equation becomes       ⅆ          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        p            ⁢              xe2x80x83            ⁢      ρ      ⁢              xe2x80x83            ⁢                        (          τ          )                *            ⁢              u        ⁡                  (                      p            ,                          t              -              px                                )                    
The inverse transform function for parabolic trajectories is as follows:
g(t,x,xcfx84,p)=xcfx84xe2x88x92t+px2
Thus, when xcfx84=txe2x88x92px2, the inverse linear slant stack Radon transformation equation becomes       ⅆ          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        p            ⁢              xe2x80x83            ⁢      ρ      ⁢              xe2x80x83            ⁢                        (          τ          )                *            ⁢              u        ⁡                  (                      p            ,                          t              -                              px                2                                              )                    
The inverse transform function for hyperbolic trajectories is as follows:
g(t,x,xcfx84,p)=xcfx84xe2x88x92{square root over (t2xe2x88x92p2x2)}
Thus, when xcfx84={square root over (t2xe2x88x92p2x2)}, the inverse hyperbolic Radon transformation equation becomes       ⅆ          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        p            ⁢              xe2x80x83            ⁢      ρ      ⁢              xe2x80x83            ⁢                        (          τ          )                *            ⁢              u        ⁡                  (                      p            ,                                                            t                  2                                -                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
The choice of which form of Radon transformation preferably is guided by the travel time trajectory at which signals of interest in the data are recorded. Common midpoint gathers, because they offer greater accuracy by measuring a single point from multiple source-receiver pairs, are preferred over other types of gathers. Primary reflection signals in a common midpoint gather generally will have hyperbolic kinematic trajectories. Thus, it would be preferable to use hyperbolic Radon transforms.
To date, however, Radon transformations based on linear slant stack and parabolic trajectories have been more commonly used. The transform function for hyperbolic trajectories contains a square root whereas those for linear slant stack and parabolic transform functions do not. Thus, the computational intensity of hyperbolic Radon transformations has in large part discouraged their use in seismic processing.
It has not been practical to accommodate the added complexity of hyperbolic Radon transformations because the computational requirements of conventional processes are already substantial. NMO correction involves a least-mean-squares (xe2x80x9cLMSxe2x80x9d) energy minimization calculation. Forward and inverse Radon transforms also are not exact inverses of each other. Accordingly, an additional LMS calculation is often used in the transformation. Those calculations in general and, especially LMS analyses, are complex and require substantial computing time, even on advanced computer systems.
Moreover, a typical seismic survey may involve hundreds of sources and receivers, thus generating tremendous quantities of raw data. The data may be subjected to thousands of different data gathers. Each gather is subjected not only to filtering processes as described above, but in all likelihood to many other enhancement processes as well. For example, data is typically processed to compensate for the diminution in amplitude as a signal travels through the earth (xe2x80x9camplitude balancingxe2x80x9d). Then, once the individual gathers have been filtered, they must be xe2x80x9cstackedxe2x80x9d, or compiled with other gathers, and subjected to further processing in order to make inferences about subsurface formations. Seismic processing by prior art Radon methods, therefore, requires substantial and costly computational resources, and there is a continuing need for more efficient and economical processing methods.
It is understood, of course, that the transformation, NMO correction, and various other mathematical processes used in seismic processing do not necessarily operate on all possible data points in the gather. Instead, and more typically, those processes may simply sample the data. For example, the transform functions in Radon transformations are based on variables t, x, xcfx84, and p. The transformations, however, will not necessarily be performed at all possible values for each variable. Instead, the data will be sampled a specified number of times, where the number of samples for a particular variable may be referred to as the index for the variable set. For example, k, l, m, and j may be defined as the number of samples in, or the index of, respectively, the time (t), offset (x), tau (xcfx84), and slowness (p) sets. The samples will be taken at specified intervals, xcex94t, xcex94x, xcex94xcfx84, and xcex94p. Those intervals are referred to as sampling variables.
The magnitude of the sampling variables, the variable domain and all other factors being equal, will determine the number of operations or samples that will be performed in the transformation and the accuracy of the transformation. As the sampling variables become finer, i.e., smaller, the number of samples increases and so does the accuracy. By making the sampling variables larger, or coarser, the number of samples is reduced, but the accuracy of the transformation is reduced as well.
In large part because of the computational intensity of those processes, however, the sampling variables used in prior art Radon processes can be relatively coarse. In NMO correction, for example, industry typically samples at xcex94t values several time greater than that at which the field data was recorded, i.e., in the range of 20-40 milliseconds, and samples at xcex94v values of from about 15 to about 120 m/sec. To the extent that the sampling variables are greater than the data acquisition sample rates, data is effectively discarded and the accuracy of the process is diminished. In the Radon transformation itself, prior art methods typically set xcex94t, xcex94xcfx84, and xcex94x at the corresponding sample rates at which the data was recorded, i.e., in the range of 2 to 4 milliseconds and 25 to 400 meters or more, respectively. xcex94p is relatively coarse, however, typically being set at values in the range of 4-24 xcexcsec/m. Also, the index j of the slowness (p) set usually is equal to the fold of the offset data, i.e., the number of source-receiver pairs in the survey, which typically ranges from about 20 to about 120.
Accordingly, prior art methods may rely on additional computations in an attempt to compensate for the loss of accuracy and resolution inherent in coarse sampling variables. For example, the LMS analysis applied in NMO correction and prior art Radon transformations are believed to obviate the need for finer sampling variables. Certain prior art processes utilize trace interpolation, thereby nominally reducing xcex94x values, but the additional xe2x80x9cdataxe2x80x9d is an approximation of what would be recorded at a particular offsets between the offsets where receivers were actually located. Relative inaccuracy in processing individual data gathers also is believed to be acceptable because the data from a large number of gathers ultimately are combined for analysis, or xe2x80x9cstackedxe2x80x9d. Nevertheless, there also is a continuing need for more accurate processing of seismic data and processes having finer resolutions.
The velocity at which primary reflection signals are traveling, what is referred to as the stacking velocity, is used in various analytical methods that are applied to seismic data. For example, it is used in determining the depth and lithology or sediment type of geological formations in the survey area. It also is used various seismic attribute analyses. A stacking velocity function may be determined by what is referred to as a semblance analysis. Which such methods have provided reasonable estimations of the stacking velocity function for a data set, given the many applications based thereon, it is desirable to define the stacking velocity function as precisely as possible.
An object of this invention, therefore, is to provide a method for processing seismic data that more effectively removes unwanted noise from meaningful reflection signals.
It also is an object to provide such methods based on hyperbolic Radon transformations.
Another object of this invention is to provide methods for removal of noise from seismic data that are comparatively simple and require relatively less computing time and resources.
Yet another object is to provide methods for more accurately defining the stacking velocity function for a data set.
It is a further object of this invention to provide such methods wherein all of the above-mentioned advantages are realized.
Those and other objects and advantages of the invention will be apparent to those skilled in the art upon reading the following detailed description and upon reference to the drawings.
The subject invention provides for methods of processing seismic data to remove unwanted noise from meaningful reflection signals indicative of subsurface formations. The methods comprise the steps of obtaining field records of seismic data detected at a number of seismic receivers in an area of interest. The seismic data comprises signal amplitude data recorded over time by each of the receivers and contains both primary reflection signals and unwanted noise events. The seismic data is assembled into common geometry gathers in an offset-time domain.
The amplitude data is then transformed from the offset-time domain to the time-slowness domain using a limited Radon transformation. That is, the Radon transformation is applied within defined slowness limits pmin and pmax, where pmin is a predetermined minimum slowness and pmax is a predetermined maximum slowness. Preferably, an offset weighting factor xn is applied to the amplitude data, wherein 0 less than n less than 1.
Such Radon transformations include the following continuous transform equation, and discrete versions thereof that approximate the continuous transform equation:       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢                        ∫                      -            ∞                    ∞                ⁢                              ⅆ            t                    ⁢                      ⅆ                          (                              x                ,                t                            )                                ⁢                      x            n                    ⁢                      δ            ⁡                          [                              f                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                          
Radon transformations based on linear slant stack, parabolic, and other non-hyperbolic kinematic travel time trajectories may be used, but those based on hyperbolic Radon kinematic trajectories are preferred. Suitable hyperbolic Radon transformations include the following continuous transform equation, and discrete versions thereof that approximate the continuous transform equation:       u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ                  xx          n                    ⁢              ⅆ                  (                      x            ,                                                            τ                  2                                +                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
Preferably, the Radon transformation is a high resolution Radon transformation. That is, the Radon transformation is performed according to an index j of the slowness set and a sampling variable xcex94p; wherein       j    =                            p          max                -                  p          min                +                  1          ⁢          μ          ⁢                      xe2x80x83                    ⁢          sec          ⁢                      /                    ⁢          m                            Δ        ⁢                  xe2x80x83                ⁢        p              ,
xcex94p is from about 0.5 to about 4.0 xcexcsec/m, pmax is a predetermined maximum slowness, and pmin is a predetermined minimum slowness.
A corrective filter is then applied to the transformed data to enhance the primary reflection signal content of the data and to eliminate unwanted noise events. Preferably, the corrective filter is a time variant, high-low filter, and it determined by reference to the velocity function derived from performing a semblance analysis or a pre-stack time migration analysis on the amplitude data.
After filtering, the enhanced signal content is inverse transformed from the time-slowness domain back to the offset-time domain using an inverse Radon transformation, and, if necessary, an inverse of the offset weighting factor pn is applied to the inverse transformed data, wherein 0 less than n less than 1. Such Radon transformations include the following continuous inverse transform equation, and discrete versions thereof that approximate the continuous inverse transform equation:       ⅆ          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        p            ⁢                        ∫                      -            ∞                    ∞                ⁢                  xe2x80x83                ⁢                              ⅆ            τ                    ⁢                      xe2x80x83                    ⁢                      p            xe2x80x3                    ⁢                                    ρ              ⁡                              (                τ                )                                      *                    ⁢                      u            ⁡                          (                              p                ,                τ                            )                                ⁢                      δ            ⁡                          [                              g                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                          
Suitable hyperbolic inverse Radon transformations include the following continuous inverse transform equation, and discrete versions thereof that approximate the continuous inverse transform equation:       ⅆ          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ                  pp          xe2x80x3                    ⁢              xe2x80x83            ⁢      ρ      ⁢              xe2x80x83            ⁢                        (          τ          )                *            ⁢              u        ⁡                  (                      p            ,                                                            t                  2                                -                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
The signal amplitude data for the enhanced signal content is thereby restored. The processed and filtered data may then be subject to further processing by which inferences about the subsurface geology of the survey area may be made.
It will be appreciated that by limiting the transformation within defined slowness limits pmin and pmax, the novel processes provide increased efficiency by reducing the amount of data transformed and increase efficacy by decimating the noise signals that fall outside those limits. In general, pmax is greater than the slowness of reflection signals from the shallowest reflective surface of interest. The lower transformation limit, pmin, generally is less than the slowness of reflection signals from the deepest reflective surface of interest. Typically, the upper and lower transformation limits pmin and pmax will be set within 20% below (for lower limits) or above (for upper limits) of the slownesses of the pertinent reflection signals, or more preferably, within 10% or 5% thereof.
It will be appreciated that primarily because they limit the Radon transformation, and avoid more complex mathematical operations required by prior art processes, the novel methods require less total computation time and resources. Thus, even though hyperbolic Radon transformations heretofore have generally been avoided because of their greater complexity relative to linear slant stack and parabolic Radon transformations, the novel processes are able to more effectively utilize hyperbolic Radon transformations and take advantage of the greater accuracy they can provide, especially when applied to common midpoint geometry gathers.
The subject invention also provides for methods for determining a stacking velocity function. The methods comprise the steps of processing and filtering the data as described above. A semblance analysis is performed on the restored data to generate a refined semblance plot, and a velocity analysis is performed on the refined semblance plot to define a refined stacking velocity function. Preferably, an offset weighting factor xn is first applied to the restored amplitude data, wherein 0 less than n less than 1.
It will be appreciated that, because the semblance and velocity analysis is performed on data that has been more accurately processed and filtered, that the resulting stacking velocity function is more accurate. This enhances the accuracy of the many different seismic processing methods that utilize stacking velocity functions. Ultimately, the increased accuracy and efficiency of the novel processes enhances the accuracy of surveying underground geological features and, therefore, the likelihood of accurately locating the presence of oil and gas deposits.